R for Geospatial Workshop
Keene Morrow
11/21/2019
Introduction to Raster Data
Raster Basics
# Load digital surface model for Harvard Forest
# 1m GEOTIF
# UTM Zone 18
# can navigatre through subfolders by hitting tab
HARV_dsmCrop_info <- GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")
HARV_dsmCrop_info## rows 1367
## columns 1697
## bands 1
## lower left origin.x 731453
## lower left origin.y 4712471
## res.x 1
## res.y 1
## ysign -1
## oblique.x 0
## oblique.y 0
## driver GTiff
## projection +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## file data/NEON-DS-Airborne-Remote-Sensing/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif
## apparent band summary:
## GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float64 TRUE -9999 1 1697
## apparent band statistics:
## Bmin Bmax Bmean Bsd
## 1 305.07 416.07 359.8531 17.83169
## Metadata:
## AREA_OR_POINT=Area# Load the GEOTIF as a raster
DSM_HARV <- raster(here::here("data", "NEON-DS-Airborne-Remote-Sensing", "NEON-DS-Airborne-Remote-Sensing", "HARV", "DSM", "HARV_dsmCrop.tif"))
DSM_HARV## class : RasterLayer
## dimensions : 1367, 1697, 2319799 (nrow, ncol, ncell)
## resolution : 1, 1 (x, y)
## extent : 731453, 733150, 4712471, 4713838 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : C:/Users/User/Documents/R/spatial-in-r_2019-11/data/NEON-DS-Airborne-Remote-Sensing/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif
## names : HARV_dsmCrop
## values : 305.07, 416.07 (min, max)# call summary of raster
summary(DSM_HARV) # using 4.31% of cells## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (4.31% of all cells)## HARV_dsmCrop
## Min. 305.32
## 1st Qu. 345.49
## Median 359.48
## 3rd Qu. 374.12
## Max. 414.43
## NA's 0.00summary(DSM_HARV, maxsamp=ncell(DSM_HARV)) # using all the cells## HARV_dsmCrop
## Min. 305.07
## 1st Qu. 345.59
## Median 359.67
## 3rd Qu. 374.28
## Max. 416.07
## NA's 0.00#DSM_HARV is currently a raster... we want a data frame to plot in ggplot
DSM_HARV_df <- as.data.frame(DSM_HARV, xy = TRUE) %>%
mutate(z = HARV_dsmCrop)
# this df has more than 2 million rows
# consider whether a raster really needs to be a data frame...
ggplot() +
geom_raster(data = DSM_HARV_df, aes(x = x, y = y, fill = z)) +
scale_fill_viridis_c() +
coord_quickmap()
plot(DSM_HARV) # makes a gross looking plot of the raster, but it's paster than ggplot for this purpose
Coordinate Reference Systems
Challenge
Use the output from the GDALinfo() function to find out what NoDataValue is used for our DSM_HARV dataset.
GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")## rows 1367
## columns 1697
## bands 1
## lower left origin.x 731453
## lower left origin.y 4712471
## res.x 1
## res.y 1
## ysign -1
## oblique.x 0
## oblique.y 0
## driver GTiff
## projection +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## file data/NEON-DS-Airborne-Remote-Sensing/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif
## apparent band summary:
## GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float64 TRUE -9999 1 1697
## apparent band statistics:
## Bmin Bmax Bmean Bsd
## 1 305.07 416.07 359.8531 17.83169
## Metadata:
## AREA_OR_POINT=AreaFound in first data.frame output from GDALinfo() NoDataValue = -9999
Histogram of raster
ggplot() +
geom_histogram(data = DSM_HARV_df, aes(HARV_dsmCrop))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Reproject Raster Data in R
Raster Projection in R
DTM_HARV <- raster(here::here("data", "NEON-DS-Airborne-Remote-Sensing", "NEON-DS-Airborne-Remote-Sensing", "HARV", "DTM", "HARV_dtmCrop.tif")) # WGS84 using UTM Zone 18
DTM_HARV## class : RasterLayer
## dimensions : 1367, 1697, 2319799 (nrow, ncol, ncell)
## resolution : 1, 1 (x, y)
## extent : 731453, 733150, 4712471, 4713838 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : C:/Users/User/Documents/R/spatial-in-r_2019-11/data/NEON-DS-Airborne-Remote-Sensing/NEON-DS-Airborne-Remote-Sensing/HARV/DTM/HARV_dtmCrop.tif
## names : HARV_dtmCrop
## values : 304.56, 389.82 (min, max)DTM_hill_HARV <- raster(here::here("data", "NEON-DS-Airborne-Remote-Sensing", "NEON-DS-Airborne-Remote-Sensing", "HARV", "DTM", "HARV_DTMhill_WGS84.tif")) # WGS84 using lat-long
crs(DTM_hill_HARV)## CRS arguments:
## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0Reprojecting the hillshade file
Reprojecting changes the data!
DTM_hill_UTMZ18N_HARV <- projectRaster(DTM_hill_HARV, crs = crs(DTM_HARV))crs(DTM_hill_UTMZ18N_HARV)## CRS arguments:
## +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0extent(DTM_hill_UTMZ18N_HARV)## class : Extent
## xmin : 731397.3
## xmax : 733205.3
## ymin : 4712403
## ymax : 4713907crs(DTM_HARV)## CRS arguments:
## +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0Math with Rasters: Canopy Height Model
Creating a canopy height model by subtracting digital terrain model from digital surface model
# double check that the projections are the same
crs(DTM_HARV)## CRS arguments:
## +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0crs(DSM_HARV)## CRS arguments:
## +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0# check that the extents are the same
extent(DTM_HARV)## class : Extent
## xmin : 731453
## xmax : 733150
## ymin : 4712471
## ymax : 4713838extent(DSM_HARV)## class : Extent
## xmin : 731453
## xmax : 733150
## ymin : 4712471
## ymax : 4713838# they are
# if they weren't, we'd use clip()
# and check them out visually
plot(DTM_HARV) # smoother and lower elevation
plot(DSM_HARV) # rougher and higher elevation #### With direct subtraction
# quick and dirty canopy height model
CHM <- DSM_HARV - DTM_HARV
plot(CHM)
#### With overlay() overlay() takes two rasters and a function telling it what to do with them
CHM_ov_HARV <- overlay(DSM_HARV, DTM_HARV, fun = function(r1, r2) {return(r1 - r2)})
plot(CHM_ov_HARV)